
proc joint_behaviour {jtype index X Y hc bc hb bb fc cv dbL dbV N H pfile} {
source Units.tcl
# $jtype 1= External node 2= Internal node
if {$jtype==2} {
set k_cr 		0.29 ; #experimental test Shear strength coefficient at cracking 
set k_pk 		0.42 ; #experimental test Shear strength coefficient at peak
set k_ult 		0.42 ; #experimental test Shear strength coefficient at ultimate  
set gamm_cr	    0.0002 ; #experimental test Shear deformation according to O'Reilly at cracking 
set gamm_pk		0.0090 ; #experimental test Shear deformation according to O'Reilly at peak
set gamm_ult	0.0200 ; #experimental test Shear deformation according to O'Reilly at ultimate  
set hyst 		[list 0.6 0.2 0.0 0.010 0.3]; #calibrated parameters from experimental test 
}
if {$jtype==1} {
set k_cr 		0.132 ; #experimental test Shear strength coefficient at cracking  
set k_pk 		0.132 ; #experimental test Shear strength coefficient at peak
set k_ult 		0.053 ; #experimental test Shear strength coefficient at ultimate  
set gamm_cr		0.0002 ; #experimental test Shear deformation according to O'Reilly at cracking
set gamm_pk		0.0132 ;  #experimental test Shear deformation according to O'Reilly at peak
set gamm_ult	0.0270 ;  #experimental test Shear deformation according to O'Reilly at ultimate 
set hyst		[list 0.6 0.2 0.0 0.0 0.3]; #calibrated parameters from experimental test 
}
# ----------------------------------------------------------------------------------------------
# ----------------------Additional parameter-----------------------------------------------------
# ----------------------------------------------------------------------------------------------
set Ac [expr $hc*$bc];
if {$bc >= $bb} {set bj [expr {min($bc,[expr $bb+0.5*$hc])}]} ;
if {$bc < $bb}  {set bj [expr {min($bb,[expr $bc+0.5*$hc])}]} ;

puts [format "Joint base: bj= %.3f" $bj]

set pt_cr [expr $k_cr*pow($fc,0.5)]; #Principle joint stress at cracking  
set pt_pk [expr $k_pk*pow($fc,0.5)]; #Principle joint stress at peak  
set pt_ult [expr $k_ult*pow($fc,0.5)]; #Principle joint stress at ultimate
set z  [expr 0.9*($hb-$cv-$dbV-$dbL/2)]; #level arm

puts [format "Principle joint stress:
pt_cr= %.3f  pt_pk= %.3f pt_ult= %.3f" $pt_cr $pt_pk $pt_ult ]
puts [format "Level arm: z= %.3f" $z]
# ----------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------
if {$jtype==2} {
# Positive Moment: Cracking peak and residual
set M_cr_p  [expr $pt_cr*$bj*$hc*(($H*1000*$z)/($H*1000-$z))*(pow(1+($N*1000/($pt_cr*$bj*$hc)),0.5))/1000000];
set M_pk_p  [expr $pt_pk*$bj*$hc*(($H*1000*$z)/($H*1000-$z))*(pow(1+($N*1000/($pt_pk*$bj*$hc)),0.5))/1000000];
set M_ult_p [expr $pt_ult*$bj*$hc*(($H*1000*$z)/($H*1000-$z))*(pow(1+($N*1000/($pt_ult*$bj*$hc)),0.5))/1000000];
# Negative Moment: Cracking peak and residual
set M_cr_n  [expr -$pt_cr*$bj*$hc*(($H*1000*$z)/($H*1000-$z))*(pow(1+($N*1000/($pt_cr*$bj*$hc)),0.5))/1000000];
set M_pk_n  [expr -$pt_pk*$bj*$hc*(($H*1000*$z)/($H*1000-$z))*(pow(1+($N*1000/($pt_pk*$bj*$hc)),0.5))/1000000];
set M_ult_n [expr -$pt_ult*$bj*$hc*(($H*1000*$z)/($H*1000-$z))*(pow(1+($N*1000/($pt_ult*$bj*$hc)),0.5))/1000000];

puts [format "Internal Joint behaviour positive:
k_cr= %.3f  k_pk= %.3f k_ult= %.3f" $k_cr $k_pk $k_ult]
puts [format "Internal Joint behaviour positive:
M_cr_p= %.3f  M_cr_p= %.3f M_ult_p= %.3f" $M_cr_p $M_pk_p $M_ult_p ]
puts [format "Internal Joint behaviour negative:
M_cr_n= %.3f  M_cr_n= %.3f M_ult_n= %.3f" $M_cr_n $M_pk_n $M_ult_n ]
}

if {$jtype==1} {
# Positive Moment: Cracking peak and residual
set a  [expr $hc/$hb*0.5]
set a2 [expr pow( $hc/$hb*0.5,2)]
set M_cr_p  [expr $pt_cr*$bj*$hc*$H*1000*$z/($H*1000-$z) * ($a + sqrt($a2 + 1 + $N*1000/($pt_cr*$bj*$hc)))/1000000];
set M_pk_p  [expr $pt_pk*$bj*$hc*$H*1000*$z/($H*1000-$z) * ($a + sqrt($a2 + 1 + $N*1000/($pt_pk*$bj*$hc)))/1000000];
set M_ult_p [expr $pt_ult*$bj*$hc*$H*1000*$z/($H*1000-$z) * ($a + sqrt($a2 + 1 + $N*1000/($pt_ult*$bj*$hc)))/1000000];
# Negative Moment: Cracking peak and residual
set M_cr_n  [expr -$pt_cr*$bj*$hc*$H*1000*$z/($H*1000-$z) * ($a + sqrt($a2 + 1 + $N*1000/($pt_cr*$bj*$hc)))/1000000];
set M_pk_n  [expr -$pt_pk*$bj*$hc*$H*1000*$z/($H*1000-$z) * ($a + sqrt($a2 + 1 + $N*1000/($pt_pk*$bj*$hc)))/1000000];
set M_ult_n [expr -$pt_ult*$bj*$hc*$H*1000*$z/($H*1000-$z) * ($a + sqrt($a2 + 1 + $N*1000/($pt_ult*$bj*$hc)))/1000000];

puts [format "Internal Joint behaviour positive:
k_cr= %.3f  k_pk= %.3f k_ult= %.3f" $k_cr $k_pk $k_ult]
puts [format "Internal Joint behaviour positive:
M_cr_p= %.3f  M_cr_p= %.3f M_ult_p= %.3f" $M_cr_p $M_pk_p $M_ult_p ]
puts [format "Internal Joint behaviour negative:
M_cr_n= %.3f  M_cr_n= %.3f M_ult_n= %.3f" $M_cr_n $M_pk_n $M_ult_n ]
}
# ----------------------------------------------------------------------------------------------
# -- Create the nodes (2 node at the same position to connect zeroLength element) --------------
# ----------------------------------------------------------------------------------------------
set $m [expr $N/9.81];
node 1${index}	$X $Y -mass $m 0.0 0.0;
node 6${index}	$X $Y -mass 0.0 0.0 0.0;
# ----------------------------------------------------------------------------------------------
# -- Create the uniaxial Hysteretic Material for joint------------------------------------------
# ----------------------------------------------------------------------------------------------
uniaxialMaterial Hysteretic 3${index} $M_cr_p $gamm_cr $M_pk_p $gamm_pk	$M_ult_p $gamm_ult -$M_cr_p -$gamm_cr -$M_pk_p -$gamm_pk -$M_ult_p -$gamm_ult [lindex $hyst 0] [lindex $hyst 1] [lindex $hyst 2] [lindex $hyst 3] [lindex $hyst 4];
# ---------------------------------------------------------------------------------------------
# -- Implement the MaxMin Limits  -------------------------------------------------------------
# ----------------------------------------------------------------------------------------------
set gamm_max 0.100; # Use a higher limit for now, needs to be updated later
uniaxialMaterial MinMax 5${index} 3${index} -min -$gamm_max -max $gamm_max
# ----------------------------------------------------------------------------------------------
# -- Create the axial material  ----------------------------------------------------------------
# ----------------------------------------------------------------------------------------------
set Ec	[expr 3320*sqrt($fc)+6900]; # Concrete Elastic Modulus 
set Kspr 	[expr 2*$Ec*$MPa*$Ac/$hb]; #axial stiffnes of the joint (2EcAc/hb)
uniaxialMaterial Elastic 1${index} 1e15; # This is for the horizontal direction 
uniaxialMaterial Elastic 2${index} $Kspr; # This is for vertical direction
# ----------------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------------
# -- Create the zeroLength elements  -----------
# ----------------------------------------------------------------------------------------------
set rigM 	1${index}; # horizontal direction 
set axM 	2${index}; # vertical direction
set flXM	5${index}; # joint flexibility due to shear
set ET	    9${index}; # element tag
element zeroLength $ET 1${index} 6${index} -mat $rigM $axM $flXM -dir 1 2 3 -doRayleigh 1
puts $pfile [format "Element 1${index} Mjx:%.2f %.2f %.2f %.2f %.2f %.2f gamma:%.2f %.2f %.2f %.2f %.2f %.2f" $M_cr_p $M_pk_p $M_ult_p $M_cr_n $M_pk_n $M_ult_n $gamm_cr $gamm_pk $gamm_ult -$gamm_cr -$gamm_pk -$gamm_ult]
}
